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Abstract.  An  investigation  into  density  and  velocity  fluctuations  in  the 
oceanic  thermocline  is  presented.  Two  kinds  of  numerical  simulation  are 
reported.  In  the  first,  an  attempt  is  made  to  capture  the  transition  from 
breaking  internal  waves  to  the  small-scale  turbulence  they  generate.  The 
model  used  for  this  is  based  on  a continual  forcing  of  a large-scale  standing 
internal- wave.  Evidence  is  presented  for  a transition  in  the  energy  spectra 
from  the  anisotropic  k~3  buoyancy  range  to  the  small-scale  fc~5/3  isotropic 
inertial  range.  Density  structures  that  form  during  wave  breaking  are  analyzed 
and  regions  of  mixing  associated  with  the  breaking  events  are  visualized.  In 
the  second  kind  of  simulation,  internal-wave  packets  are  followed  as  they 
propagate  through  the  thermocline.  It  is  found  that  the  breaking  of  crests 
within  the  packet  can  lead  to  overturning  events  on  the  scale  observed  in  the 
ocean,  and  the  subsequent  turbulence  can  form  a continuous  wake. 


Introduction 

In  recent  observations  of  fluctuations  in  the  oceanic 
thermocline,  Alford  and  Pinkel  (1999,  2000)  found  many 
overturns  with  vertical  scale  of  about  2 m and  these 
were  highly  correlated  with  the  presence  of  energetic 
waves  with  vertical  wavelengths  on  the  order  of  10  m. 
Large  scale  fluctuations,  say  10  m and  above  in  ver- 
tical scale,  can  be  described  reasonably  well  as  inter- 
nal waves.  For  much  smaller  scales,  say  1 m and  be- 
low, the  flow  is  probably  better  described  in  terms  of 
nearly  isotropic  turbulence.  Intermediate  between  the 
large-scale  wave  dynamics  and  the  small-scale  turbu- 
lence is  a transition  regime  in  which  there  is  a com- 
petition between  waves  and  turbulence.  It  is  this 
intermediate  range,  often  called  the  buoyancy  range, 
that  contains  the  overturning  activity  observed  by  Al- 


ford and  Pinkel.  Since  the  observations  are  essen- 
tially one-dimensional  in  space,  a direct  numerical  sim- 
ulation which  could  faithfully  describe  events  in  this 
range  would  help  toward  understanding  the  full  three- 
dimensional  flow  structures  behind  the  observations. 
The  first  part  of  our  investigation  will  focus  on  the 
production  of  overturns  by  an  idealized  internal  wave 
forcing. 

The  second  part  of  our  investigation  concerns  the 
propagation  of  internal- wave  packets  through  the  ther- 
mocline. Alford  and  Pinkel  (1999,  2000)  observed  co- 
herent regions  of  strong  oscillatory  vertical  strain  rate 
that  travel  vertically  through  100  m or  more  of  the  ther- 
mocline. These  propagating  structures  had  an  internal 
wave  structure  with  vertical  wavelength  of  about  10  m, 
and  the  entire  coherent  region  could  be  described  by 
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an  envelope  of  about  50  m in  vertical  extent.  These 
coherent  structures  were  strongly  correlated  with  2 m 
overturns.  Given  the  complicated  nature  of  the  flow, 
with  strong  components  of  large-scale  advection,  it  was 
difficult  to  know  precisely  what  kind  of  structures  to 
associate  with  these  coherent  localized  oscillations.  Al- 
ford and  Pinkel  suggested  that  these  were  internal- wave 
packets.  Here  we  hope  to  validate  that  identification 
to  a certain  extent  by  showing  how  an  idealized  packed 
would  propagate  through  a simplified  model  of  the  ther- 
mocline,  and  by  showing  that  the  observed  overturn- 
ing scales  could  be  consistent  with  internal-wave  packet 
propagation. 

Given  current  computer  resources,  it  would  not  be 
practical  to  simulate  all  of  the  scales  that  are  relevant 
to  the  breaking  of  waves  in  the  buoyancy  range.  The 
forcing  of  the  waves  is  thought  to  result  from  a com- 
plicated interaction  of  many  internal  waves  with  scales 
ranging  in  the  vertical  from  10’s  of  meters  to  hundreds 
of  meters  and  more,  and  typical  horizontal  internal  wave 
scales  can  be  even  much  greater  than  these.  On  the 
small-scale  end  of  the  spectrum,  the  observed  break- 
ing is  occurring  on  vertical  scales  of  the  order  of  1 m, 
and  these  breaking  events  produce  turbulence  that  ex- 
tends down  to  a viscous  cutoff  on  the  order  of  1 cm. 
Thus  direct  numerical  simulation  of  the  entire  range  of 
scales  is  still  impractical.  We  will  use  a combination 
of  LES  modeling  and  an  artificial  model  of  the  large 
scale  forcing  in  order  to  reduce  the  spectral  range  that 
we  will  need  to  cover.  To  attack  our  first  problem  of 
investigating  how  waves  at  the  short-scale  end  of  the 
Garrett-Munk  (1975)  spectrum  go  unstable  and  break 
in  the  buoyancy  range,  we  have  used  an  artificial  forcing 
with  length-scales  fixed  at  20  m in  the  vertical  and  20 
m in  the  horizontal  to  represent  the  effect  of  all  larger 
scales.  At  the  small-scale  end  of  the  simulation,  we  have 
introduced  an  eddy  viscosity  with  a cutoff  at  the  16  cm 
level  in  both  horizontal  and  vertical  directions.  Thus 
our  model  does  some  violence  to  the  true  physics  at 
the  large  and  small-scale  ends  of  the  simulated  range. 
However,  the  hope  is  that  it  will  do  justice  to  the  evo- 
lution in  the  buoyancy  range.  This  model  does  prove 
capable  of  capturing  the  transition  from  the  buoyancy 
to  the  inertial  range.  For  the  problem  of  wave-packet 
propagation,  which  is  perhaps  the  source  of  the  order 
10  m scale  variability  most  correlated  with  overturning, 
we  needed  to  expand  our  domain  size  in  order  to  allow 
for  the  propagation  and  evolution  of  the  packet. 

Forced  20-meter  wave 

Our  first  goal  is  to  determine  to  what  extent  our 
simulations  can  capture  the  transition  from  the  buoy- 
ancy range  to  the  inertial  range  in  the  energy  spec- 


tra. Constructing  a theory  of  this  transition  is  compli- 
cated because  of  the  anisotropy  of  the  buoyancy  range. 
To  make  progress,  some  theoretical  formulations  have 
represented  the  entire  spectrum  as  depending  only  on 
wavenumber  k.  The  model  for  the  kinetic  energy  spec- 
trum in  the  buoyancy  range  is  then 


E(k)  = aN2k 


2u~3 


(1) 


where  a is  an  empirical  constant  and  N is  the  Brunt- 
Vaisala  frequency,  which  measures  the  strength  of  the 
stratification.  The  Brunt-Vaisala  frequency  is  defined 

by 


N2^_9_dp 
Po  dz 


(2) 


where  g is  the  acceleration  of  gravity,  p is  the  back- 
ground density  profile,  assumed  stable  (i.e.  dp/dz  < 0), 
and  po  is  the  volume  average  of  p.  From  the  observed 
spectra  of  vertical  shear,  the  constant  a is  determined 
to  be  about  0.47,  but  it  will  be  more  convenient  for  us 
to  consider  the  two  components  of  the  horizontal  veloc- 
ity ( u , v)  separately,  and,  assuming  horizontal  isotropy 
in  the  observations,  this  would  suggest  a « 0.2  for  the 
spectrum  of  either  component  (cf.  Gibson  1986,  Gargett 
et  al.  1981).  The  inertial  range  kinetic  energy  spectrum 
is  given  by 

E(k)  = CKe2/3k~5/3  (3) 

where  e is  the  turbulent  dissipation  rate  of  total  kinetic 
energy  and  Ck  is  the  empirical  Kolmogorov  constant. 
A reasonable  value  to  assume  for  Ck  is  1.5  (cf.  Lesieur, 
1997).  For  the  energy  of  one  component  of  the  velocity 
field,  there  would  simply  be  a prefactor  of  1/3  multi- 
plying this  isotropic  spectrum.  The  Ozmidov  (or  buoy- 
ancy) wavenumber  is  then  estimated  by  simply  match- 
ing these  two  spectra  at  wavenumber  kb-  The  result  is, 
up  to  an  order  one  multiplicative  constant, 


kb  = y/N*/e, 


(4) 


(cf.  Holloway,  1981;  Gibson  1986). 

The  model  for  the  potential  energy  spectrum  in  the 
buoyancy  range  is  similar  to  that  for  the  kinetic  en- 
ergy spectrum.  The  empirical  constant  a for  the  tem- 
perature spectrum  is  found  to  have  value  of  about  0.2 
(cf.  Gibson  1986,  Gregg  1977).  The  spectral  model  for 
the  inertial  range  of  density  fluctuations  is  the  Corrsin- 
Obukhov  spectrum,  which  involves  the  decay  rate  of 
density  fluctuations  as  well  as  e.  For  our  purposes,  we 
prefer  to  write  the  spectrum  directly  in  terms  of  the 
turbulent  decay  rate  of  potential  energy,  which  we  shall 
write  as  epe.  Then  the  Corrsin- Obukhov  spectrum  for 
the  potential  energy  takes  the  following  form: 

PE(k)  = C0epee-1/3fc-5/3, 


(5) 
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where  C0  is  the  Corrsin  constant. 

For  all  of  the  simulations  presented  here,  we  have 
used  the  Boussinesq  approximation.  The  evolution 
equations  may  be  written  as 

~ + u-  Vu+  — Vp'  - — g = i/(V2)V2u,  (6) 

at  po  po 

V • u = 0.  (7) 

l + + (8) 

where  i/(-)  and  «(■)  are  considered  functions  of  the 
Laplacian  operator  and  are  used  to  represent  eddy  pa- 
rameterizations  in  general  (cf.  Herring  and  Metais 
1992)  and  g = — gz.  We  have  neglected  the  effect  of  ro- 
tation, which  should  not  play  a major  role  at  the  small 
scales  with  which  we  are  concerned.  The  total  density 
is  given  by 

p = p(z)  + p'(x,y,z,t),  (9) 

where  p'(x,  y,  z,t)  is  the  deviation  from  the  horizona- 
tally  averaged  density  p(z).  po  is  the  average  of  p(z) 
over  2.  The  pressure  p'  is  the  deviation  from  the  back- 
ground mean  pressure,  p'  can  be  determined  in  terms  of 
u by  taking  the  divergence  of  (6)  under  the  assumption 
that  the  velocity  field  is  divergenceless. 

We  simulate  these  dynamical  equations  with  a spec- 
tral code  with  triply  periodic  boundary  conditions.  As  a 
sub-grid  scale  parameterization,  we  have  used  the  large- 
eddy  simulation  model  of  Lesieur  and  Rogallo  (1989). 
This  eddy  viscosity  i /t(k)  is  approximately  constant 
throughout  the  buoyancy  range  and  the  large-scale  end 
of  the  inertial  range,  but  increases  rapidly  with  k in  the 
vicinity  of  the  spectral  cutoff  kc.  Due  to  the  spectral 
shape  of  the  eddy  viscosity,  this  model  is  sometimes 
called  the  cusp  model.  It  seems  reasonable  in  modeling 
the  buoyancy  range  to  use  such  a model  since  it  does 
not  completely  neglect  the  effects  of  unresolved  eddies 
on  the  buoyancy  range,  but,  at  the  same  time,  it  puts 
the  strongest  eddy  viscosity  in  the  inertial  range  near 
the  cutoff. 

We  should  emphasize  the  point  that  the  size  of  the 
eddy  viscosity  depends  on  the  amount  of  energy  at  the 
cutoff  scale.  If  the  resolution  of  the  simulation  of  a given 
physical  flow  is  increased,  that  is  if  kc  is  increased,  then 
the  eddy  viscosity  will  be  correspondingly  smaller.  The 
total  viscosity  used  in  the  simulations  is  the  sum  of 
the  eddy  viscosity  and  the  constant  molecular  viscosity 
vmoi.  Thus  the  z/(V2)  in  equation  (6)  in  the  spectral 
simulation  is  taken  as  the  total  viscosity: 

v{k)  = Vmoi  +Vt{k)  (10) 


the  closure  model  for  stratified  turbulence.  For  simplic- 
ity, we  have  just  taken  the  turbulent  Prandtl  number 
Prt(k)  to  be  a fixed  constant  independent  of  k in  our 
simulations.  We  determined  this  constant  by  examin- 
ing the  evolution  of  the  potential  energy  spectrum  for 
decaying  stratified  turbulence  that  is  initially  highly  ex- 
cited at  all  scales.  More  specifically,  we  started  with  an 
initial  spectrum  in  which  the  GM  spectrum  was  contin- 
ued to  scales  below  10  m as  in  the  decay  simulations 
of  Siegel  and  Domaradzki  (1994).  With  Prt  = 0.55 
our  simulations  of  decaying  turbulence  produced  spec- 
tra with  the  high  wavenumbers  obeying  the  fc~5/3  law 
for  both  velocity  and  density  fluctuations. 

Next  we  turn  to  the  question  of  the  forcing.  The 
large-scale  flows  that  actually  drive  the  buoyancy  range 
are  predominantly  the  waves  of  the  Garrett-Munk  range. 
The  full  range  where  internal  wave  dynamics  dominates 
includes  scales  of  kilometers  in  the  horizontal  and  hun- 
dreds of  meters  in  the  vertical.  Because  of  lack  of  reso- 
lution, we  cannot  provide  a full  representation  of  the 
effects  of  all  large-scale  internal  wave  forcing  on  the 
buoyancy  range.  In  our  model,  of  necessity,  we  perform 
a drastic  reduction  in  modeling  the  forcing;  we  replace 
the  driving  of  all  of  the  GM  waves  by  a linear  stand- 
ing wave  at  one  wavelength.  Bouruet-Aubertot  et  al. 
(1995,  1996),  in  two-dimensional  simulations  of  a strat- 
ified turbulence,  excited  a standing  wave  of  the  type 
we  use,  but  they  allowed  this  wave  to  decay,  whereas 
we  maintain  its  amplitude  at  the  same  level  throughout 
the  simulation. 

To  give  the  form  of  the  forcing  used,  let  us  first  in- 
troduce nondimensional  units.  We  will  take  all  lengths 
to  be  scaled  by  27 t/L,  where  L is  the  length  of  one  side 
of  our  computational  domain.  Time  will  be  scaled  by 
1/N.  Then  the  frequency  of  linear  internal  waves  is 
given  by 

»-±x-  (11) 

where  kh  = + fc2  is  the  horizontal  wavenumber. 

One  particular  linear  standing  wave  is 


* £ 
u = ( u,v,w ) — A-~(0,  sin  y sin  z,  cosy  cos  z)  sin  —j=, 

(12) 

— = A cos  u cos  z cos  — =,  (13) 

Po  y/2 

where  A is  an  arbitrary  amplitude  and  g*  is  the  nondi- 
mensional gravity.  Note  that  the  dimensional  period  of 
this  wave,  which  is  the  forcing  period,  is  given  by 


r~~  27T 

Tr  = ^ N' 


(14) 


The  choice  of  turbulent  diffusion  depends  on  the 
choice  of  values  for  various  parameters  that  enter  into 


To  give  some  idea  of  the  structure  of  this  standing 
wave,  we  show  in  Figure  1 a contour  plot  of  the  density 
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Figure  1.  Contour  plot  of  the  density  field  in  a vertical 
y-z  cross  section  through  the  center  of  the  domain.  The 
width  and  height  of  the  cross  section  are  each  20  m, 
which  corresponds  to  2-tt  in  nondimensional  units.  The 
instant  shown  corresponds  to  maximum  displacement 
of  the  isopycnals  for  the  forced  standing  wave. 


field  in  a vertical  y — z cross  section.  Note  that  the  den- 
sity field  in  this  standing  wave  has  no  variation  along 
the  x direction.  In  this  figure,  we  see  an  instantaneous 
representation  of  the  iso-density  surfaces.  When  t/  \/2  is 
an  odd  multiple  of  7r/2,  these  isopycnals  will  all  be  flat. 
The  degree  to  which  they  deviate  from  that  at  other 
times  depends  on  the  value  of  A as  well  as  t.  The  instant 
of  time  represented  here  is  such  that  t/V 2 is  an  integer 
multiple  of  7r  and,  hence,  one  of  maximum  distortion  of 
the  density  contours.  Note  that  the  density  field  in  (13) 
has  two  nodal  planes,  represented  by  two  thick  contour 
lines  in  the  figure,  at  z — ±tt/2  (nondimensional).  Dur- 
ing the  forcing  cycle,  these  planes  remain  flat  and  fixed 
in  position.  The  fluid  above  and  below  these  planes  ver- 
tically approaches  and  retreats  from  them  depending  on 
the  phase  in  time  and  the  y position  considered.  Thus 
the  points  on  the  nodal  planes  at  y = 0 and  y = ±tt 
are  the  centers  of  regions  of  oscillating  high  strain  rate. 
On  the  other  hand  the  points  where  the  isopycnals  are 
steepest,  that  is  where  y = ±7r/2,  ±37t/2  and  z — 0,  ±n 
the  magnitude  of  the  shear  dv/dz  is  highest.  Thus  one 
advantage  of  the  standing  wave  forcing  is  that  the  points 
of  highest  shear  and  highest  strain  rate  remain  fixed  in 
space  making  it  easier  to  differentiate  the  kinds  of  over- 
turning events  associated  with  these  extremes.  This  will 
be  convenient  for  making  comparisons  with  Alford  and 


Pinkets  (2000)  analysis.  In  particular,  they  noted  that 
there  were  regions  of  high  shear  where  classical  shear 
instabilities  often,  but  not  always,  resulted  in  overturn- 
ing. Even  more  more  interestingly,  there  were  regions  of 
high  Ri  in  which  overturns  were  also  observed.  In  more 
than  half  of  these  cases  Ri  was  even  greater  than  2 sug- 
gesting that  the  typical  shear  instability  (Ri  < 1/4)  is 
unlikely.  Many  of  these  overturns  were  in  regions  of  high 
vertical  strain  rate.  With  our  standing  wave  forcing,  the 
regions  of  high  shear  and  high  strain  rate  are  separate 
and  each  occur  in  the  same  location  during  each  forcing 
oscillation.  This  helps  simplify  the  analysis. 

In  two-dimensional  numerical  studies  of  Bouruet- 
Aubertot  et  al.  (1996),  the  standing  wave  becomes  un- 
stable and  generates  turbulence.  This  would  also  hap- 
pen in  our  three-dimensional  simulation,  but  the  tur- 
bulence would  be  highly  constrained  since  there  is  yet 
no  source  of  x-variation  in  our  flow.  To  break  the  two- 
dimensional  symmetry  of  the  flow,  while  maintaining 
the  basic  structure  of  the  large  scale,  we  add  a weak 
component  of  forcing  with  ^-variation.  We  have  tried 
this  in  various  ways:  adding  a random  initial  pertur- 
bation at  all  scales,  randomly  forcing  the  modes  with 
k = 1 at  each  time  step,  adding  another  large-scale 
standing  wave,  adding  a propagating  wave,  and  so  on. 
The  results  are  similar  to  each  other  if  the  perturbations 
are  sufficiently  weak.  For  the  simulations  discussed  be- 
low, we  have  added  to  the  primary  forcing  wave  only  a 
small  amplitude  standing  wave  of  the  same  spatial  scale. 
Specifically,  we  added  the  following  perturbation: 


u = j4'-p(cos(a:  + z),0,  - cos(x  + z ))  sin  —j=,  (15) 

v2  v2 


t 

Po 


A'  cos(x  + z)  cos  — . 

V2 


(16) 


Thus  in  the  simulations  discussed  below  the  forcing  oc- 
curs only  at  k — \fl.  The  coefficient  A!  was  taken  to 
be  .4/20,  and,  hence,  the  energy  in  the  perturbation  is 
only  1/400  that  of  the  primary  forcing  wave. 

We  performed  a series  of  experiments  in  which  the 
size  of  the  computational  domain  and  the  amplitude  of 
the  forcing  were  varied.  The  initial  studies  were  at  res- 
olution 643  and  showed  that  for  sufficiently  large  am- 
plitudes A for  which  the  forcing  wave  itself  was  over- 
turning, a A:-5/3  spectrum  extending  over  most  of  the 
spectral  range  could  be  established.  For  weaker  forc- 
ing, a steeper  spectrum  approximating  k~3  was  found 
( Camevale  and  Briscolini , 1999).  For  intermediate  am- 
plitude forcings,  we  were  able  to  observe,  at  least  inter- 
mittently, cases  which  do  appear  to  exhibit  the  transi- 
tion from  the  buoyancy  range  to  the  smaller  scale  in- 
ertial range.  Weak  and  strong  forcings  are  measured 
relative  to  shear  amplitudes  typical  in  the  thermocline. 
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The  best  results  were  obtained  with  a forcing  ampli- 
tude that  could  actually  be  considered  representative 
of  wave  amplitudes  in  the  thermocline.  Specifically,  the 
forcing  amplitude  that  we  refer  to  as  intermediate,  is  for 
a value  of  A in  equations  (12)  and  (13)  such  that  the 
maximum  shear  during  a cycle  of  the  forcing  is  equiva- 
lent to  the  rms  shear  of  the  GM  spectrum  at  the  scale  of 
our  computational  domain.  The  rms  shear  is  calculated 
by  integrating  the  shear  of  the  GM  spectrum  from  the 
kilometer  scale  down  to  the  scale  of  interest  (cf.,  Gregg, 
1989).  Our  best  results  tended  to  be  for  cases  in  which 
the  vertical  wavelength  of  the  forcing  was  20  m.  For 
N = 3 cph,  the  net  rms  shear  from  the  GM  spectrum 
for  this  scale  is  Sgm  (20  m)«  3 x 10-3s-1  (cf.,  Gregg, 
1989).  Taking  this  value  to  determine  the  amplitude 
of  our  forcing,  we  obtain  a standing  wave  in  which  the 
largest  deviation  of  the  density  isosurfaces  are  as  illus- 
trated in  Figure  1.  Thus  we  have  a standing  wave  that 
does  not  itself  overturn  during  the  forcing  cycle,  and, 
in  addition,  the  Richardson  number  of  the  forcing  wave, 
defined  by 
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does  not  drop  below  3.125.  Therefore,  the  forcing  wave 
itself  is  convectively  stable  and  not  subject  to  shear  in- 
stability. This  kind  of  forcing  is  consistent  with  the  pic- 
ture that  the  GM  waves  themselves  are  not  convectively 
or  shear  unstable,  but  through  wave- wave  interactions 
will  produce  smaller  scale  waves  that  are  unstable  by 
these  criteria.  Choosing  a stronger  forcing  wave  that 
is  itself  convectively  or  shear  unstable  would  miss  the 
important  cascade  process  that  produces  the  unstable 
waves  of  the  buoyancy  range,  but  would  rather  produce 
turbulence  directly  resulting  in  an  inertial  range  (cf. 
Camevale  and  Briscolini , 1999). 

For  all  of  the  simulations  discussed  below,  we  used 
a resolution  1283  and  a computational  cube  of  20  m on 
a side.  Our  isotropic  spectral  cutoff  is  at  wavenumber 
60,  and  the  smallest  resolved  wavelength  is  about  33 
cm  (with  grid  spacing  20  m/128  « 16  cm).  The  forc- 
ing amplitude  was  fixed  so  that  the  max  shear  in  the 
forced  wave  is  Sgm  (20  m),  and  the  Vaisala  frequency 
was  taken  to  be  3 cycles  per  hour,  which  is  a typical 
oceanic  value. 

A long  simulation  was  performed  with  realistic  values 
for  the  molecular  viscosity  and  diffusivity.  The  kine- 
matic viscosity  was  set  to  vmoi  = 0.01  cm 2/s  and  the 
molecular  Prandtl  number  at  Prmoi  = 7 (cf.  Gargett, 
1985).  We  can  calculate  a Reynolds  number  for  the 
oceanic  flow  for  vertical  motions  on  the  20  m scale  by 
using  the  rms  shear.  Thus  we  can  write 


Re  = Sgm{L)L2 /vmoi-  (18) 


For  L — 20  m,  this  Reynolds  number  would  be  ap- 
proximately 105.  By  including  the  molecular  viscosity, 
the  simulation  is  an  attempt  to  represent  flow  with  this 
Reynolds  number.  We  will  see  that  there  is  not  much 
difference  with  results  obtained  by  neglecting  the  molec- 
ular components  of  viscosity  and  diffusivity.  That  is  to 
say  that  over  the  range  of  scales  simulated  (20  m to  33 
cm)  the  difference  between  infinite  Reynolds  number 
flow  and  that  for  Re  = 10s  is  small. 

We  can  think  of  our  standing  wave  forcing  as  the 
linear  susperposition  of  a set  of  propogating  internal 
waves.  To  be  precise,  the  combination  of  the  two  stand- 
ing waves  given  in  (12)  and  (15)  consists  of  12  propagat- 
ing plane  waves.  These  wave  interact  nonlinearly  pro- 
ducins  smaller-scales  that  eventually  fill  out  the  entire 
spectrum.  The  early  evolution  is  essentially  just  that  of 
the  nearly  two-dimensional  standing  wave.  During  this 
time  there  are  only  sinusoidal  waves  on  the  most  dis- 
turbed isosurface,  but  these  waves  then  fold  over  form- 
ing elongated  overturns.  These  regions  are  convectively 
unstable  and  break.  At  this  point  the  three  dimension- 
ality of  the  flow  becomes  evident. 

After  about  five  cycles  of  the  forcing,  the  large-scale 
wave  breaks  repeatedly,  however,  not  necessarily  during 
each  forcing  period.  The  wave  breaking  on  the  most 
disturbed  isopycnal  occurs  roughly  symmetrically  with 
large-scale  overturning  occurring  nearly  at  the  same  val- 
ues of  y and  z each  time  and  along  lines  of  constant  x, 
respecting  in  the  large  scales  the  symmetry  of  the  main 
part  of  the  forcing.  However,  no  two  breaking  events 
with  the  subsequent  evolution  during  the  forcing  cy- 
cle are  the  same.  In  Figure  2,  for  one  such  cycle,  we 
show  eight  instantaneous  images  of  this  isosurface  us- 
ing a perspective  three-dimensional  plot.  The  frames 
are  ordered  temporally  from  left  to  right  and  top  to 
bottom.  The  first  frame  in  the  upper  left  hand  cor- 
ner corresponds  to  t = 11.392V,  and  the  interval  be- 
tween frames  is  At  = Tp/7.  Thus  the  first  and  last 
frames  correspond  to  the  same  phase  of  the  forcing. 
The  first  frame  captures  the  moment  when  breaking 
is  just  beginning.  Let  us  say  that  the  first  four  frames 
represent  the  breaking  event,  and  the  last  four  the  af- 
termath. We  see  that  during  the  breaking  event,  heavy 
fluid  spills  over  lighter  fluid,  crashing  down  with  the  cre- 
ation of  small-scale  structures  all  along  the  fines  of  the 
two  breaking  regions.  Similar  behavior  is  in  laboratory 
experiments  with  standing-wave  forcing  ( Taylor  1992, 
McEwan  1983a).  Afterwards,  the  region  of  the  small- 
scale  turbulent  structures  spreads,  eventually  ’contam- 
inating’ the  entire  isosurface.  If  we  compare  the  final 
frame  with  the  first  frame,  we  see  that  the  final  surface 
is  much  rougher,  filled  with  small-scale  structures  ev- 
erywhere, and  that  there  is  no  larger  scale  folding-over 
of  the  surface  as  there  was  in  the  first  frame.  In  the 
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F igure  2.  The  evolution  of  the  p — po  isopycnal  during  one  cycle  of  the  forcing.  The  frames  are  ordered  by  time 
from  left  to  right  and  top  to  bottom.  The  first  corresponds  to  t = 11.39 Tp,  and  the  interval  between  frames  is 
At  = TF/7. 
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later  evolution,  the  wave  will  break  again,  but  only  af- 
ter a refractory  period,  in  this  case  of  about  two  forcing 
cycles. 

Next  we  will  consider  the  energy  spectra  for  the  flow 
at  the  same  times  as  those  illustrated  in  Figure  2.  Since 
the  energy  is  highly  anisotropic  at  scales  larger  than 
those  in  the  inertial  range,  plotting  the  total  energy  as 
a function  of  the  isotropic  wavenumber  tends  to  obscure 
the  transition  between  small  and  large  scales.  To  most 
clearly  display  the  transition,  we  have  found  it  useful  to 
consider  the  spectrum  Ev(k)  of  v,  the  y component  of 
velocity,  which  is  the  horizontal  component  that  is  di- 
rectly affected  by  the  forcing.  Along  with  the  spectra, 
we  have  also  drawn  lines  corresponding  to  the  inertial 
range  spectrum  (l/3)Cife2/3fc~5/3  and  buoyancy  range 
spectrum  0.2N2k~3.  For  each  frame,  e is  taken  as  the 
total  kinetic  energy  dissipation  rate  at  that  time.  We 
have  included  a factor  of  (1/3)  which  is  appropriate 
for  a single  component  in  the  isotropic  inertial  range. 
For  the  Kolmogorov  constant,  a value  of  1.5  was  used 
in  each  case.  We  should  emphasize  that  no  attempt 
is  made  here  to  fit  the  data,  but  the  coefficient  is  just 
taken  as  this  standard  value  a priori.  For  the  buoyancy 
range  spectrum,  we  have  used  the  coefficient  a = 0.2 
in  all  cases.  In  each  frame  shown  in  Figure  3,  we  see  a 
fairly  good  match  at  wavenumbers  greater  than  about 
20  (that  is  for  scales  below  about  1 m)  to  the  Kol- 
mogorov inertial  range  spectrum.  The  main  deviation 
is  at  wavenumbers  near  k = 60,  the  cutoff  wavenumber, 
and  this  is  to  be  expected  from  previous  experience  with 
the  cusp  model  (cf.  Lesieur  and  Rogallo,  1989).  The 
spectrum  below  wavenumber  20  is  naturally  far  more 
irregular  than  that  above  due  to  the  much  smaller  num- 
ber of  modes  in  the  lower  spectral  bands.  If  we  neglect 
the  first  few  wavenumbers,  then  there  is  some  evidence 
here  for  a steeper  spectral  range  for  wavenumbers  below 
about  k = 20,  that  is  for  scales  larger  than  about  1 m, 
at  least  in  the  frames  that  correspond  to  times  during 
the  breaking  of  the  wave  (first  four  panels).  In  the  af- 
termath of  breaking,  the  spectra  tend  to  be  somewhat 
flatter  (the  last  four  panels).  The  best  representative 
of  the  transition  between  buoyancy  and  inertial  range 
is  found  in  panel  (c),  which  corresponds  to  a time  when 
the  enstrophy  is  near  a local  maximum.  Here  the  buoy- 
ancy range  spectrum  makes  a reasonably  good  fit  in 
the  range  of  scales  from  about  4 m down  to  about  1 m. 
From  the  forcing  scale  (20  m vertical)  to  about  the  5 m 
scale,  there  is  a dip  in  the  energy  that  has  also  been  seen 
in  the  spectra  from  similar  two-dimensional  simulations 
of  the  decay  of  a standing  wave  ( Bouruet-Aubertot  et  al. 
1996). 

In  this  experiment  it  appears  that  the  expected  spec- 
tral signature  of  a transition  between  a buoyancy  range 
at  large  scale  and  the  inertial  range  at  small  scale  oc- 


curs only  for  periods  during  which  there  is  active  break- 
ing. Indeed,  it  appears  that  wave-wave  interactions  re- 
peatedly build  up  energy  in  the  buoyancy  range  until 
a fc-3  spectrum  is  achieved.  At  that  point  significant 
breaking  occurs  and  energy  drains  from  the  buoyancy 
range.  Let  us  focus  on  the  breaking  event.  In  Fig- 
ure 4a  is  the  image  of  the  p — p0  isosurface  at  the 
time  identified  as  the  best  for  illustrating  the  spectral 
transition  from  the  buoyancy  to  the  inertial  range.  It 
shows  the  curling  over  and  spilling  down  or  plunging  of 
the  heavier  fluid  over  lighter,  while  Figure  4b  suggests 
mixing  by  the  appearance  of  many  small-scale  struc- 
tures along  the  two  parallel  lines  of  the  breaking  wave. 
The  corresponding  spectra  for  all  three  components  of 
kinetic  and  for  the  potential  energy  are  shown  in  Fig- 
ure 5.  First  we  notice  that  although  the  spectra  are 
highly  anisotropic  from  the  forcing  scale  (20  m)  down 
to  about  the  1 m scale,  there  is  an  approximate  ‘re- 
turn’ to  isotropy  for  the  smaller  scales.  This  is  particu- 
larly evident  in  the  kinetic  energy  spectra  for  t =11.82 
Tp  (panel  c).  In  panels  (a)  and  (c),  we  have  made  an 
attempt  to  draw  the  best  fit  inertial  range  spectra  to 
determine  the  appropriate  Kolmogorov  constants  (Ck) 
that  fit  these  data.  We  did  this  for  the  Ev  (fc)  spectra, 
obtaining  the  best  fit  ‘by  eye’  from  enlarged  portions 
of  the  small  scale  spectra.  The  result  that  was  used  to 
draw  the  inertial  range  model  spectra  in  panels  (a)  and 
(c)  is  (Ck)  = 1-4.  In  panels  (b)  and  (d),  the  poten- 
tial energy  spectra  are  drawn.  In  these  panels  the  small 
scales  were  fit  to  the  Corrsin- Obukhov  spectrum  to  de- 
termine the  appropriate  Corrsin  constant.  In  panels  (b) 
and  (d)  the  Corrsin  constants  used  to  draw  the  model 
Corrsin-Obukhov  spectrum  were  C0  = 0.83  and  0.8  re- 
spectively. In  all  panels  the  model  buoyancy  range 
spectrum  drawn  is  0.2 N2k~3.  Thus  the  Kolmogorov 
constant  found  here  is  somewhat  smaller  than  the  em- 
pirical values  of  1.5  and  the  Corrsin  constant  is  some- 
what larger  than  the  empirical  value  of  0.67.  Never- 
theless, the  values  are  remarkably  close  to  the  empiri- 
cal values,  given  that  the  spectral  width  of  the  inertial 
range  here  only  covers  wavelengths  from  about  1 me- 
ter to  about  33  cm.  Also  the  near  collapse  of  the  three 
kinetic  energy  spectra  for  small  scales  is  encouraging. 
Thus  it  seems  that  the  subgrid  scale  model  is  work- 
ing well  at  small  scales  and  that  the  dynamics  of  the 
transition  from  anisotropic  buoyancy  to  the  isotropic 
inertial  range  is  acting  as  imagined  in  theoretical  mod- 
els. Finally,  we  should  note  that  the  value  of  e from  the 
simulations  is  about  one  third  of  the  value  observed  by 
Alford  and  Pinkel  (2000)  associated  with  values  of  N = 
3 cph.  This  appears  quite  reasonable  given  the  level 
of  modeling  we  have  had  to  employ  for  the  forcing  and 
subgrid  scale  vortices. 

Besides  the  kinetic  and  potential  energy  spectra,  we 
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Figure  4.  A breaking  event  visualized  on  the  p = po  isopycnal.  These  are  enlargements  of  the  images  shown  in 
the  composite  Figure  2 in  panels  3 and  4,  corresponding  to  times  (a)  11.68  and  (b)  11.82  T p (one  seventh  of  a 
forcing  period  apart). 
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Figure  5.  (a)  Kinetic  energy  spectra  for  all  three  components  of  the  velocity  at  time  t = 11.687’ p.  The  thick  long 
dashed,  solid  and  short  dashed  lines  correspond  to  the  energy  spectra  for  the  u,v,  and  w components  respectively. 
The  thin  solid  lines  correspond  to  the  Kolmogorov  spectrum  (1/3)C p-e2^3k~5^3  with  Ck  = 1.4  and  the  saturation 
spectrum  0.2 N2k~3.  (b)  Potential  energy  spectrum  at  time  t = 11.687V.  The  thick  solid  line  corresponds  to  the 
potential  energy  spectrum.  The  thin  solid  lines  correspond  to  the  Corrsin- Obukhov  spectrum  C'0ep<.e_1/3fc"'5/3  with 
C0  = 0.83  and  the  saturation  spectrum  0.2N2k~3.  (c)  As  in  (a)  but  for  t=11.82  Tp  and  Ck  = 1.4.  (d)  As  in  (b) 
but  for  t=l  1.82  Tp  and  C0  = 0.80. 
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can  also  find  predictions  for  the  buoyancy  flux  spectrum 
in  both  the  theory  of  Lumley-Shur  (cf.  Lumley  1964, 
1967,  Phillips  1967,  Weinstock  1985)  and  the  theory 
of  Holloway  (1983,  1986).  The  modal  spectrum  of  the 
buoyancy  flux  can  be  written  as 

-gU  < wl,p'k  > / po-  (19) 

If  this  quantity  is  positive,  then  for  wavevector  k there 
is  conversion  of  potential  energy  to  kinetic  energy,  and 
vice  versa  if  it  is  negative. 

The  prediction  of  the  Lumley-Shur  theory  for  the 
buoyancy  flux  spectrum  in  the  buoyancy  and  inertial 
ranges  is 

BF(k)  = -2 (l  + D{kb/k)A^  V2  (kb/k)7/3  (20) 

where  kb  is  as  defined  in  (4)  and  D is  a constant.  Lumley 
(1964)  assumed  the  buoyancy  flux  to  be  negative  and, 
hence,  D to  be  positive.  In  displaying  his  final  result, 
Lumley  incorporated  D into  his  definition  for  kb,  but 
we  will  leave  it  explicit.  Lumley’s  prediction  of  negative 
buoyancy  flux  through  the  buoyancy  and  inertial  ranges 
is  just  the  opposite  of  what  we  have  found  numerically 
for  our  wave-forced  problem.  All  of  the  ingredients  for 
an  alternative  prediction  of  the  buoyancy  flux  are  given 
in  Holloway  (1983),  and  based  on  this  we  have  derived 
the  same  prediction  as  given  in  (20,  but  with  the  sign 
of  D clearly  arbitrary  (for  details  see  Camevale  et  al., 
2001). 

In  Figure  6a,  we  plot  the  buoyancy  flux  spectrum 
from  our  simulation  as  a function  of  k.  This  is  a time 
averaged  spectrum,  where  we  have  averaged  over  a pe- 
riod of  67V,  with  time  increment  of  0.17V.  The  time 
averaging  is  necessary  to  remove  temporal  fluctuations 
in  the  large-scales.  Note  that  the  buoyancy  flux  spec- 
trum is  negative  for  large  scales  (1  < k < 3),  and  pos- 
itive for  smaller  scales.  This  implies  a transformation 
of  kinetic  to  potential  energy  at  large  scales  (closest  to 
the  forcing  scale  k = \JT)  and  a transfer  of  potential 
to  kinetic  energy  at  all  smaller  scales.  Since  our  ob- 
served buoyancy  flux  spectrum  is  positive  through  both 
the  buoyancy  and  inertial  ranges,  it  can  be  compared  to 
the  theoretical  prediction  given  by  (20)  only  by  choos- 
ing a negative  value  for  D.  To  define  the  constant  D, 
we  note  that  the  wavenumber  where  the  buoyancy  flux 
vanishes  is  determined  by  D.  Here  we  shall  choose  D 
so  that  the  zero  value  occurs  at  k = 3.5  (corresponding 
in  our  simulation  to  a wavelength  of  5.7  m)  since  our 
buoyancy  flux  was  found  to  vanish  between  k = 3 and 
k = 4.  The  theory  will  apply  only  above  this  wavenum- 
ber, and  we  can  think  of  this  as  the  lower  limit  on  the 
buoyancy  range,  or  the  upper  wavenumber  of  the  Gar- 
rett Munk  spectrum  in  the  schematic  shown  in  our  in- 
troduction. To  compute  kb,  given  by  (4),  we  use  the 


Figure  6.  (a)  Graph  of  the  buoyancy  flux  spectrum 
averaged  over  6 periods  of  the  forcing,  with  10  samples 
per  period,  (b)  as  in  (a)  but  only  the  positive  portion  of 
the  spectrum  plotted  in  log-log  format  to  compare  with 
the  theoretical  spectrum  of  equation  (20)  with  negative 
coefficient  D.  The  result  from  the  simulation  is  repre- 
sented by  the  thick  line,  while  the  theoretical  spectrum, 
based  on  eo  = e where  the  overbar  represents  time  av- 
eraging, is  drawn  as  a thin  line.  All  graphs  in  (a)  and 
(b)  are  normalized  by  e/Kb. 
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time  averaged  dissipation  rate  e.  Thus  all  parameters 
in  the  theory  are  determined  by  the  wavenumber  where 
the  buoyancy  flux  vanishes  and  the  values  of  N and  l 
(in  this  case  kb  = (Ar3/e)1/2  * 34.1).  The  resulting  the- 
oretical buoyancy  flux  spectrum  is  compared  in  Figure 
6b  to  the  results  from  our  simulation.  For  wavenum- 
bers in  the  buoyancy  range,  the  match  between  theory 
and  simulations  is  reasonably  good.  For  the  theoreti- 
cal curve,  the  decay  with  k is  approximately  k~ 7 ^ for 
all  k above  about  10.  The  simulation  data  follow  the 
theoretical  curve  fairly  well  up  to  about  wavenumber 
20,  where  the  simulation  spectrum  begins  to  deviate 
from  k~ 7/3,  and  is  clearly  much  shallower  than  this  for 
k > 30.  This  shallowness  of  the  simulation  spectrum 
for  k greater  than  about  30  is  probably  an  indication 
that  the  buoyancy  flux  is  not  captured  properly  by  the 
SGS  model  near  the  high  wavenumber  cutoff.  The  cusp 
model  viscosity  grows  rapidly  with  k for  wavenumbers 
above  about  k = 30  and  is  largest  at  kmax ■ This  is 
just  the  range  where  our  buoyancy  flux  spectrum  be- 
comes very  shallow.  It  is  very  possible  that  the  artificial 
damping  of  the  high  k modes  that  the  model  performs 
to  mimic  transfer  of  energy  beyond  kmax  does  not  allow 
for  the  proper  treatment  of  the  buoyancy  flux  in  that 
region.  But  this  is  not  unexpected  for  such  a subgrid- 
scale  model. 

Positive  buoyancy  flux  for  small  scales  has  also  been 
found  in  other  simulations.  In  direct  numerical  simula- 
tions (i.e.  simulations  without  subgrid  scale  modeling) 
of  forced  stratified  turbulence  in  both  two  and  three- 
dimensions,  Holloway  (1988),  and  Ramsden  and  Hol- 
loway (1992)  showed  that  the  buoyancy  flux  was  nega- 
tive only  at  large  scales  and  positive  at  small  scales. 
These  results  were  interpreted  as  meaning  negative 
buoyancy  flux  for  k < kb  (i.e.  in  the  buoyancy  range) 
and  positive  buoyancy  flux  for  higher  k.  However,  the 
forcing  used  in  their  simulations  was  spectrally  fairly 
broad,  and  it  would  not  be  inconsistent  with  their  re- 
sults to  say  that  the  buoyancy  flux  was  negative  at  the 
strongly  forced  modes  and  positive  for  smaller  scales  as 
in  our  findings.  Additionally,  we  have  repeated  our  nu- 
merical simulations  with  a finite  difference  code  using 
a Smagorinsky  eddy  viscosity,  a very  independent  test, 
and  also  found  positive  buoyancy  flux  through  the  buoy- 
ancy and  inertial  ranges.  In  their  finite  difference  LES 
study  of  shear  driven  stratified  turbulence,  Kaltenbachet 
et  al.  (1994)  also  found  positive  buoyancy  flux  at  small 
scales  and  negative  at  large  scales,  although  we  must 
note  that  the  region  of  negative  buoyancy  flux  in  their 
simulations  is  spectrally  very  broad  compared  to  ours. 
In  two-dimensional  flow  simulations  of  the  decay  of  a 
standing  wave  of  just  the  type  that  we  use  for  forcing 
our  flow,  Bouruet-Aubertot  et  al.  (1996)  found  that  the 
buoyancy  flux  was  positive  through  most  of  the  range 


that  they  identified  with  the  buoyancy  range,  and  also 
that  the  flux  followed  a fc~7/3  spectral  law  in  a run  with 
grid  resolution  256 2 and  a slightly  steeper  law  at  reso- 
lution 5122  (note  that  those  simulations  did  not  include 
an  inertial  range). 

There  was  some  discussion  at  the  ’Aha  Huliko’a 
meeting  about  the  possibility  that  the  sign  of  the  buoy- 
ancy flux  found  in  these  simulations  is  affected  by  the 
type  of  forcing  used.  Our  forcing  inputs  both  potential 
and  kinetic  energy,  while  it  may  be  more  suitable  to  con- 
sider a source  of  kinetic  energy  alone  in  these  problems. 
If  there  were  no  explicit  external  forcing  of  the  density 
evolution,  then  the  net  buoyancy  flux  (averaged  over 
time)  would  have  to  be  negative  to  balance  the  drain 
of  potential  energy  due  to  diffusion.  This  however  does 
not  mean  that  the  buoyancy  flux  would  have  to  be  nega- 
tive for  all  scales.  Further  simulations  would  be  helpful 
to  define  how  the  buoyancy  flux  spectrum  varies  as  the 
mix  of  kinetic  and  potential  energy  sources  are  changes. 

Structures  in  regions  of  high  strain  rate 

The  main  structure  of  interest  in  the  buoyancy  range 
evident  in  the  density  isosurfaces  presented  in  the  last 
section  is  the  overturn  produced  by  the  curling  over  of 
the  isosurface  in  a manner  familiar  from  surface  wave 
breaking.  The  overturning  region  shown  in  the  breaking 
wave  illustrated  in  Figure  2c  has  a vertical  scale  of  about 
2 meters.  This  is  similar  in  size  to  overturns  found 
in  oceanographic  measurements  in  the  buoyancy  range. 
Alford  and  Pinkel  (2000)  made  an  inventory  of  more 
than  2200  overturns.  They  found  a median  Thorpe 
scale,  a measure  of  the  vertical  extent  of  the  overturn,  of 
1.88  m.  Note  that  this  is  not  greatly  different  from  the 
scale  suggested  by  the  transition  point  in  the  spectra 
shown  in  Figure  5,  where  the  associated  length  scale  is 
about  1.2  m.  Since  the  observational  data  are  primarily 
one-dimensional  in  space,  it  is  difficult  to  form  a three- 
dimensional  image  of  those  overturns.  The  ability  to 
perform  three-dimensional  analysis  of  such  structures 
is  one  of  the  benefits  of  numerical  simulation. 

Examining  the  full  density  field  more  thoroughly,  we 
also  find  interesting  structures  of  a rather  different  na- 
ture than  those  associated  with  strong  vertical  shear. 
These  can  be  represented  well  by  the  deformations  of 
the  density  surfaces  that  are  the  flat  nodal  surfaces  of 
the  forcing  wave.  We  shall  just  refer  to  these  surfaces  as 
the  nodal  surfaces  even  when  perturbed  and  deformed 
by  eddies.  The  most  basic  motion  of  the  fluid  in  the 
nodal  surfaces  is  alternately  toward  and  away  from  the 
centers  of  high  strain;  however,  the  combination  of  the 
large-scale  background  straining  motion  and  small-scale 
eddies  produces  localized  deformations  of  the  nodal  sur- 
face that  can  result  in  overturning  and  mixing  in  a man- 
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Figure  7.  The  evolution  of  three  density  isosurfaces  ( g*(p  — Po)/po  = 0,  n/2,  jr)  showing  the  evolution  of  ‘spouts’ 
from  a ‘nodal  surface’  and  their  subsequent  collapse  with  considerable  broadening  and  mixing.  Times  represented 
are  t = 12.1,  12.4,  12.5,  12.7,  12.8,  and  13.1  TF. 
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ner  different  from  the  overturns  discussed  in  the  last 
section.  By  plotting  simultaneously  three  density  iso- 
surfaces (one  ‘nodal  surface’  and  the  most  strongly  per- 
turbed isosurfaces  above  and  below  it)  for  a sequence  of 
times  during  the  forcing  cycle,  we  can  get  some  under- 
standing of  the  nature,  formation,  evolution,  and  fate 
of  these  structures.  In  Figure  7,  the  sequence  proceeds 
from  left  to  right,  top  to  bottom.  We  have  shifted  our 
view  of  the  computational  domain  by  an  amount  in  the 
vertical  sufficient  to  center  the  upper  ‘nodal  surface’  in 
the  image.  Above  and  below  the  ‘nodal  surface,’  the 
most  strongly  displaced  isosurfaces  are  shown.  The  up- 
per and  lower  isosurfaces  move  vertically  but  always 
in  opposite  directions  at  any  horizontal  location.  The 
combined  effect  of  the  motion  of  these  surfaces  above 
and  below  the  ‘nodal  surface’  produces  vertical  ‘dila- 
tion’ and  ’compression’  centered  on  the  ‘nodal  surface’ 
without  producing  large-scale  sinusoidal  displacement 
of  that  surface.  In  panels  (a),  (b)  and  (c),  the  up- 
per/lower surface  is  moving  upward/downward  in  the 
middle  of  the  domain  (i.e.  at  y — 0,  where  y is  the  hor- 
izontal coordinate),  and  oppositely  at  the  left  and  right 
ends  of  the  domain.  This  is  associated  with  the  verti- 
cal straining  of  the  nodal  surface  in  the  middle  and  at 
the  left  and  right  ends  of  the  domain.  In  panel  (b)  the 
isosurface  ‘erupts’  with  elements  moving  up  and  down 
along  a midline  pointing  into  the  plane.  The  eruption 
reaches  its  maximum  extension  when  the  upper  and 
lower  surfaces  stop  their  motion,  and  reverses  direction 
around  the  time  of  panel  (c).  At  y — 0 on  the  ‘nodal 
surface’  this  is  a time  of  maximum  vertical  dilational 
strain  but  zero  strain  rate  (where  dwfdz  is  the  vertical 
strain  rate).  The  structures  formed  by  these  eruptions 
represent  localized  intrusions  of  heavy  fluid  into  light 
fluid  and  vice  versa.  We  shall  refer  to  them  as  ‘spouts.’ 
As  the  upper  and  lower  isosurfaces  move  back  toward 
the  ‘nodal  surface,’  the  sense  of  straining  motion  is  re- 
versed and  the  spouts  that  were  formed  are  flattened. 
This  causes  a spreading  out  of  these  structures,  which 
in  some  cases  results  in  tossing  elements  of  the  spouts 
to  the  right  and  left  of  the  midline.  This  leads  to  the 
kind  of  pattern  seen  in  panel  (e)  which  is  in  part  an 
elongated  horizontal  structure  as  opposed  to  the  elon- 
gated vertical  structures  originally  produced  during  the 
vertically  dilational  phase  of  the  large-scale  straining 
motion.  The  final  panel  (f)  shows  the  isosurface  a short 
time  after  the  upper  and  lower  surfaces  have  again  re- 
versed their  direction  of  vertical  motion.  This  is  a phase 
of  the  motion  near  to  that  of  the  initial  panel  (a),  but 
now  there  is  a mixed  patch  of  fluid  at  the  mid  section 
(y=0)  of  the  ‘nodal  surface.’ 

To  summarize,  we  can  say  that  the  ‘spouts’  origi- 
nate from  small-scale  deviations  of  the  nodal  surface 
created  by  turbulent  flow  at  the  nodal  surface.  Once 


perturbations  pull  structures  from  the  nodal  planes  ver- 
tically, these  elements  are  subject  to  advection  due  to 
the  large-scale  straining  motion  of  the  forcing  wave.  At 
times  and  positions  where  the  straining  is  highly  di- 
lational in  the  vertical,  these  deviations  from  the  flat 
plane  elongate  vertically  and  narrow  horizontally,  form- 
ing ‘spouts.’  Then,  during  the  vertically  compressional 
and  horizontally  dilational  phase  of  the  forcing,  the 
spout  is  elongated  horizontally  creating  regions  of  con- 
vert ively  unstable  overturned  fluid.  Note  that  if  the 
large-scale  forcing  were  the  only  field  acting  on  the 
spout,  than  the  growth  of  the  spout  would  simply  have 
been  reversed  when  the  sense  of  the  straining  motion 
was  reversed.  Thus  the  presence  of  the  eddy  field  must 
play  an  important  role  in  this  irreversible  process.  The 
distortions  of  the  spout  by  the  eddy  field  are  enhanced 
during  the  horizontally  dilational  phase  of  the  evolu- 
tion. 

Internal-wave  packets 

The  observations  of  Alford  and  Pinkel  (2000)  show 
vertically  propagating  structures  at  depths  from  150  to 
350  m which  they  suggest  may  be  internal  wave  pack- 
ets. These  structures  have  vertical  extent  of  about  50  m 
with  internal  vertical  wavelengths  of  about  12  m and  are 
associated  with  overturning  events  with  vertical  scales 
of  about  2 m.  Recent  theoretical  analysis  by  Thorpe 
(1999)  provides  a criterion  for  determining  whether  the 
small-scale  turbulence  generated  by  the  overturns  in  a 
packet  will  be  left  behind  in  just  small  patches  or  in  con- 
tinuous ‘scars’  much  longer  than  the  size  of  the  packet. 
Stimulated  by  these  developments,  we  have  embarked 
on  a numerical  investigation  of  internal  wave  packets. 

Assuming  a constant  background  Brunt-Vaisala  fre- 
quency N and  ignoring  the  effects  of  the  earth’s  ro- 
tation, the  intrinsic  dimensional  frequency  for  internal 
waves  is 

(21) 

The  observed  frequency  for  one  of  the  wavepackets  in 
the  Alford  and  Pinkel  (2000)  data  is  4 cph.  This  is 
higher  than  the  ambient  N « 3 cph.  Since  amax  — N,  it 
is  assumed  that  the  observed  frequency  for  this  packet 
is  the  sum  of  the  intrinsic  frequency  plus  a Doppler 
shift.  To  predict  this  shift,  it  is  necessary  to  know  the 
wavelength  of  the  packet,  the  magnitude  of  the  ambient 
current  and  its  direction  relative  to  the  packet  propa- 
gation direction.  Alford  and  Pinkel  (2000)  suggest  that 
the  intrinsic  frequency  for  their  packet  with  observed 
frequency  of  4 cph  is  near  0.14  cph  which  leads  one  to  a 
wavelength  of  180  m.  This  suggests  that  the  horizontal 
wavelengths  in  both  directions  are  much  larger  than  the 
vertical  wavelength.  For  our  numerical  modeling,  this 
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represents  a difficulty.  We  are  reluctant  to  introduce 
anisotropic  grids  for  fear  of  the  distortions  that  might 
result,  especially  when  applying  simple  sub-grid  scale 
models.  Thus,  in  this  preliminary  work,  we  decided  to 
consider  only  the  case  in  which  horizontal  and  vertical 
wavelengths  were  equal.  The  corresponding  intrinsic 
frequency  would  then  be  about  2 cph  which  would  still 
be  consistent  with  the  observed  packet,  just  requiring 
less  of  a Doppler  shift  to  match  the  observed  frequency. 

As  for  the  amplitude  of  the  observed  packets,  this 
can  be  given  in  terms  of  the  peak  magnitude  of  the  ob- 
served strain  rate  dw/dz.  The  maximum  value  of  verti- 
cal strain  rate  in  the  Alford  and  Pinkel  (2000)  observa- 
tions is  approximately  N,  and  in  the  case  of  the  particu- 
lar packet  discussed  above,  it  seems  that  the  maximum 
is  about  0.381V. 

In  what  follows,  we  will  examine  the  evolution  of 
a particular  wave  packet  with  two-dimensional  simula- 
tions. In  an  attempt  to  reproduce  the  kind  of  behavior 
evident  in  the  observations,  we  used  simulations  in  a 
domain  of  200  m in  both  width  and  depth.  We  used 
a packet  with  non-dimensional  wavenumbers  of  12  in 
both  directions,  corresponding  to  vertical  and  horizon- 
tal wavelengths  of  (200  m)/12  as  17  m.  Our  2D  simula- 
tions had  an  effective  resolution  corresponding  to  a cut- 
off wavelength  of  as  0.8  m.  To  follow  this  phenomenon 
in  DNS  with  all  relevant  scales  well  resolved  would  re- 
quire resolution  from  200  m down  to  a few  cm,  which 
is  somewhat  impractical.  Since  the  subgrid  scale  model 
used  in  the  3D  simulations  is  not  appropriate  in  2D, 
we  had  recourse  to  hyperviscosity  (with  the  Laplacian 
taken  to  the  eighth  power).  The  simulations  illustrated 
here  are  from  a spectral  code  dealiased  with  the  3/2 
rule  ( Orszag  1971).  Although  there  are  768  wavevec- 
tors  used  in  each  direction,  after  application  of  the  3/2 
rule  this  leaves  only  512  active  modes  in  each  direction. 

Linear  dispersion  of  packets 

The  linearized  version  of  the  Boussinesq  evolution 
equations  can  be  used  to  obtain  a model  of  the  inter- 
nal wave  packet.  The  vorticity  and  density  of  a plane 
internal  wave  can  be  written  dimensionally  as 

(ojx,u>y,u>z,p')  = Aekexpi(k  • r - at),  (22) 

where  A is  an  arbitrary  amplitude  and  e is  the  eigen- 
vector 

®k  = (jjkky/Nkh,  gkkx j Nkh,0,  po) • (23) 

Taking  a linear  superposition  of  such  waves  distributed 
continuously  in  wavevector  space  and  centered  on  a par- 
ticular wavevector,  say  ko,  would  produce  an  internal 


wave  packet.  For  example, 

(«,  p’)  = KeJ  G(k  - k0)ekei(k  r-fft)d3fc,  (24) 

with 

G{ p)  = A exp  (~a2p2x  - b2p2y  - c2p2z) , (25) 

where  a,  b,  and  c are  length  scales,  represents  a propa- 
gating ellipsoidal  packet.  A slight  generalization  based 
on  simple  coordinate  rotations  will  also  permit  an  arbi- 
trary choice  for  the  orientation  of  the  ellipsoidal  enve- 
lope relative  to  the  crests  internal  to  the  packet.  Within 
the  envelope,  the  vorticity  and  density  fields  will  have  a 
phase  velocity  in  the  direction  of  ko  and  group  velocity 

cg  = Vkcrk,  (26) 

which  is  perpendicular  to  the  phase  velocity. 

By  varying  the  dimensions  a,  b and  c,  we  can  change 
the  shape  of  the  packet  as  needed.  A likely  candidate 
for  the  packets  whose  effects  are  observed  in  Alford  and 
Pinkel’s  (2000)  data  would  suggest  that  at  least  one  of 
these  lengthscales  is  very  large.  For  the  present  calcu- 
lations we  take  a to  be  infinite.  Then  we  chose  b and 
c and  the  orientation  of  the  system  to  be  such  that  the 
envelope  is  an  ellipse  with  major  axis  aligned  along  the 
direction  of  propagation.  Other  choices  may  also  be  of 
interest,  but  that  will  be  explored  in  future  work.  With 
the  ellipse  as  chosen,  the  phase  velocity  is  directed  along 
the  short  axis  and  the  group  velocity  along  the  long 
axis.  In  a numerical  simulation,  the  packet  can  only 
be  approximated,  with  the  integral  replaced  by  a dis- 
crete sum  of  wavevectors.  By  using  (24)  and  (25)  with 
t = 0,  we  are  able  to  construct  the  initial  condition  for 
a packet  that  is  both  reasonably  confined  in  space  and 
well  resolved  internally. 

The  first  issue  that  we  need  to  address  is  the  disper- 
sive spreading  of  the  wave  packet.  Simple  arguments 
suggest  that  the  physical  extent  of  the  wave  packet  will 
grow  as  Acgt  in  the  direction  of  the  group  velocity, 
where  A cg  represents  the  spread  in  group  velocities  cal- 
culated for  the  individual  wavevectors  that  contribute 
significantly  to  the  wave  packet.  We  can  make  some 
crude  dimensional  estimates  for  the  rate  of  dispersion 
by  setting  cg  ~ JV/fco  and  A cg  ~ (N/kfyAko,  where 
Ako  measures  the  spread  of  wavenumbers  in  the  packet. 
If  we  call  Axq,  the  initial  length  of  the  wavepacket,  then 
the  change  in  the  size  of  the  packet  can  be  crudely  taken 
as 

Ax  - Axo  = A cgt.  (27) 

The  packet  would  then  double  in  size  by  a time  £<*  ~ 
Axo/A cg,  and  the  distance  that  the  packet  can  travel 
before  doubling  is 


x/Ax0  k0/ Ako- 


(28) 
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Figure  8.  Contours  of  the  magnitude  of  the  pertur- 
bation density  \p'  / Po\  from  a simulation  of  the  linear 
propagation  of  a wavepacket.  The  domain  size  is  200 
m on  each  side.  The  vertical  axis  is  depth.  The  only 
contour  level  drawn  is  that  at  0.5  of  the  maximum  field 
value.  The  time  sequence  of  the  panels  is  (a)  t=0,  (b) 
t=3  hr,  (c)  t=6  hr,  and  (d)  11  hr. 


For  the  packet  used  in  the  simulations  this  predicts  a 
doubling  after  propagation  of  about  200  m. 

In  Figure  8,  we  show  the  evolution  of  the  density 
perturbation  field  during  the  propagation  of  our  packet 
following  purely  linear  dynamics.  In  each  panel,  only 
the  contour  level  corresponding  to  0.5|p'/pol  is  drawn. 
Positive  and  negative  values  have  not  been  indicated, 
but  clearly  the  sign  of  p'  will  alternate  from  one  wave 
crest  to  the  next.  We  see  the  packet  propagates  along 
the  diagonal.  This  is  in  agreement  with  the  fact  that  the 
wavevector  is  k = (12, 12)  and  that  the  group  velocity 
is  perpendicular  to  this.  It  is  less  obvious  from  the 
few  panels  that  we  can  include  here  that  the  phase  of 
the  waves  within  the  packet  advances  in  the  direction 
of  k.  The  average  speed  of  the  packet  in  propagating 
from  one  corner  of  the  domain  to  the  opposite  comer 
is  correctly  given  by  |cs|.  Furthermore,  we  see  that  the 
width  and  length  of  the  packet  grow  to  a little  more 
than  double  their  original  values  in  the  time  it  takes  to 
cross  from  one  corner  of  the  domain  to  the  other,  and 
this  is  correctly  predicted  by  the  formula  (28).  During 
the  period  of  evolution  illustrated,  the  peak  amplitude 
of  the  packet  decays  to  25%  of  its  initial  value. 

Although  the  amplitude  of  the  packet  can  be  changed 


Figure  9.  Contours  of  p/po  from  a simulation  of  the 
linear  propagation  of  a wavepacket.  The  domain  size  is 
200  m on  each  side.  The  vertical  axis  is  depth.  The  time 
sequence  of  the  panels  is  the  same  as  in  Figure  8.  The 
contour  increment  is  such  that  the  vertical  separation 
between  unperturbed  isopycnals  is  8 m. 

arbitrarily  in  this  purely  linear  simulation,  we  may  sim- 
ply assign  an  amplitude  to  see  the  effect  of  such  a packet 
on  the  full  density  field.  This  is  done  in  Figure  9.  The 
amplitude  used  represents  fluctuations  in  dw/dz  about 
five  times  the  maximum  actually  observed  in  the  Alford 
and  Pinkel  (2000)  data.  Nevertheless,  we  have  used  this 
packet  with  exaggerated  amplitude  to  more  clearly  il- 
lustrate the  nature  of  the  linear  propagation.  In  such 
a strong  packet,  there  are  regions  of  strong  overturn- 
ing, which,  if  the  packet  is  not  propagating  too  rapidly, 
would  develop  convective  instability  under  the  full  non- 
linear dynamics. 

Nonlinear  propagation  of  packets 

Having  determined  that  our  packet  propagates  cor- 
rectly under  linear  dynamics,  we  then  investigated  its 
evolution  with  the  complete  Boussinesq  equations.  The 
amplitude  of  the  observed  packet  discussed  in  the  intro- 
duction is  such  that  the  maximum  value  of  the  strain 
rate  dw/dz  is  about  0.38-/V.  With  the  packet  amplitude 
set  to  match  this  value  as  its  maximum  dw/dz,  we  per- 
formed the  simulation  illustrated  by  contour  plots  of 
p' / p in  Figure  10.  This  figure  should  be  compared  to 
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Figure  10.  Contours  of  the  magnitude  of  the  pertur- 
bation density  |p'/Po|  from  a simulation  of  the  nonlinear 


propagation  of  a wavepacket  with  max  dw/dz  w 0.38 N.  Figure  11.  Contours  of  p/p0  from  the  same  simulation 

The  domain  size  is  200  m on  each  side.  The  vertical  axis  as  represented  in  Figure  10.  The  domain  size  is  200 

is  depth.  The  only  contour  level  drawn  is  that  at  0.5  m on  each  side.  The  vertical  axis  is  depth.  The  two 

of  the  maximum  field  value.  The  time  sequence  of  the  times  illustrated  correspond  to  the  first  and  last  times 

panels  is  the  same  as  in  Figure  8.  of  Figure  10.  The  contour  increment  is  such  that  the 

vertical  separation  between  unperturbed  isopycnals  is  8 
m. 


the  corresponding  figure  for  linear  evolution,  Figure  8. 
The  times  represented  are  the  same  in  each  figure.  By 
the  time  of  panel  (b)  a clear  asymmetry  in  the  form 
of  the  packet  has  developed  in  the  nonlinear  case  and 
there  is  some  clear  distortion  of  the  packet  in  the  fi- 
nal panel.  Nevertheless,  the  overall  evolution  of  this 
nonlinear  packet  is  not  very  different  from  the  linear 
case.  This  packet  is  so  weak  that  the  initial  condition 
is  not  overturning  anywhere  and  the  Richardson  num- 
ber is  above  1 everywhere.  Thus,  the  classical  criteria 
for  convective  instability  and  shear  instability  are  not 
satisfied  in  this  packet.  This  continues  to  be  the  case 
throughout  the  simulation  in  spite  of  small-scale  gen- 
eration by  nonlinear  wave-wave  interactions.  An  idea 
of  how  weak  this  packet  is  can  be  obtained  graphically 
from  the  plots  of  the  density  contours  as  illustrated  in 
Figure  11. 

The  next  case  that  we  will  treat  is  one  for  which  the 
amplitude  of  the  packet  is  just  above  the  threshold  for 
overturning.  The  amplitude  of  this  packet  in  terms  of  its 
maximum  strain  rate  is  dw/dz  = 0.76 N.  In  Figure  12, 


we  display  the  contour  plots  for  the  perturbation  den- 
sity at  the  same  times  as  in  the  previous  figures.  We  see 
that  there  is  some  early  production  of  small  scales  that 
are  evident  in  the  wake  of  the  packet.  By  t = 6 hr  the 
packet  itself  has  become  badly  distorted,  and  by  t = 11 
hr,  it  has  degenerated  into  small-scale  structures,  al- 
though these  still  retain  to  some  extent  an  organization 
and  alignment  related  to  the  original  structure  of  the 
packet.  To  better  illustrate  the  decay  of  this  packet, 
we  display  contour  plots  of  the  full  density  field  from 
t = 2.5  hr  to  t = 4.8  hr  in  Figure  13.  Each  frame  is  an 
enlarged  image  centered  on  the  wave  packet,  showing 
only  a portion  of  the  domain  (a  square  of  size  200/3 
m on  a side).  In  panel  (a)  we  see  an  early  stage  in 
which  the  wave  is  overturning  at  points,  but  there  has 
not  yet  been  any  strong  production  of  energy  in  scales 
smaller  than  2 m (note  that  the  spacing  between  the 
unperturbed  isopycnals  is  2 m).  There  are  four  rela- 
tively strong  crests  evident  in  panel  (a).  These  crests 
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Figure  12.  Contours  of  the  magnitude  of  the  pertur- 
bation density  \p'/po\  from  a simulation  of  the  nonlinear 
propagation  of  a wavepacket  with  max  dw/dz  « 0.76 N. 
The  domain  size  is  200  m on  each  side.  The  vertical  axis 
is  depth.  The  only  contour  level  drawn  is  that  at  0.5 
of  the  maximum  field  value.  The  time  sequence  of  the 
panels  is  the  same  as  in  Figure  8. 

are  advancing  from  bottom-left  to  top-right  in  these  fig- 
ures. The  weakest  crest  (bottom-left)  is  just  entering 
the  packet  in  panel  (a).  In  the  linear  evolution  as  each 
crest  passes  through  the  packet  from  bottom- left  to  top- 
right,  its  amplitude  first  increases  and  then  decreases. 
As  envisioned  by  Thorpe  (1999),  the  crests  amplify  as 
they  move  toward  the  center  of  the  packet  and  break 
leaving  small-scale  perturbations  behind  that  link  up 
with  the  ‘debris’  produced  by  the  passage  of  previous 
crests.  The  period  of  the  sequence  of  panels  shown  here 
is  long  enough  for  the  weak  crest  on  the  lower-left  side  of 
the  packet  in  panel  (a)  to  move  completely  through  the 
packet,  finally  becoming  the  weak  crest  on  the  upper- 
right  side.  In  the  case  we  have  simulated  here,  the  crests 
do  produce  overlapping  zones  of  small-scale  perturba- 
tions that  form  a somewhat  continuous  scar,  a possibil- 
ity suggested  by  Thorpe  (1999).  One  should  note,  how- 
ever, that  during  the  period  when  a particular  crest  is 
actually  breaking,  the  overturning  and  small-scale  pro- 
duction is  not  uniform  along  the  length  of  the  crest, 
as  assumed  in  Thorpe’s  idealized  model,  but  rather  ap- 
pears in  spots  along  the  crest  (see  panels  (c)  and  (d)). 
Also  the  breaking  and  subsequent  scar  formation  does 


not  continue  indefinitely.  The  strength  of  the  packet  is 
both  dispersed  and  dissipated,  so  that  by  t = 210AT"1 
the  process  of  scar  formation  has  ceased. 

We  have  also  performed  3D  simulations  of  the  prop- 
agation of  these  wave  packets.  The  general  evolution 
exhibited  in  the  2D  simulations  is  also  found  in  3D, 
although  in  3D  we  did  not  have  sufficient  resolution  ad- 
equately capture  the  2 m overturns.  Further  details  can 
be  found  in  Camevale  and  Orlandi  (2000). 

Theory  of  packet  evolution 

We  have  addressed  here  questions  about  the  longevity 
of  wave  groups.  A theory  for  the  evolution  of  the  packet 
envelope  has  been  developed  by  Shrira  (1981)  via  mul- 
tiscale analysis  in  both  space  and  time.  His  result  for 
the  evolution  of  the  amplitude  A of  the  packet  in  the 
two-dimensional  case  studied  above  is 

iAr  + {cTkykyAyy  + 2<TkvkzAyz  + (?kzkz  Azz} 

2 

— kykyAyyy  + 3<Tjtvfcv  kZ  AyyZ  + 3crfcv  kzkxAyzz 
+vkzkzkzAzzz}  = -*7A(AA;  - A* As),  (29) 

where  ky  and  kz  are  components  of  the  central  wavector 
of  the  packet,  a is  the  intrinsic  frequency  corresponding 
to  the  central  wavevector,  s is  the  coordinate  in  the 
direction  of  propagation  of  the  packet,  d/dr  — d/dt  + 
Cgd/ds,  cg  is  the  magnitude  of  the  group  velocity,  and 

7 = k3/(akykz).  (30) 

The  typical  equation  that  arises  for  the  evolution  of  a 
wave-packet  envelope  is  the  cubic  Schroedinger  equa- 
tion which  is  significantly  different  from  (29).  It  turns 
out  that  the  term  corresponding  to  the  nonlinearity  in 
the  cubic  Schroedinger  vanishes  identically  here  due  to 
the  fact  that  a plane  wave  cannot  interact  with  itself. 
Thus  Shrira  had  to  go  to  third  order  in  the  multiple 
scale  analysis  to  obtain  the  first  contributions  of  the 
nonlinearity  to  the  evolution.  This  still  involves  a cubic 
term  for  the  nonlinearity,  but  now  not  the  simple  A|A|2 
of  the  cubic  Schroedinger  equation,  and  the  presence  of 
the  third  order  spatial  derivatives  from  the  linear  terms 
further  complicates  matters.  Notice  that  the  coefficients 
depend  in  magnitude  and  sign  on  the  orientation  of  the 
central  wavevector  of  the  packet.  Thus  we  can  antici- 
pate interesting  results  as  this  wavevector  is  varied.  In 
addition,  there  are  nonlocal  nonlinear  terms  that  arise  if 
the  flow  is  three  dimensional  that  further  greatly  com- 
plicate the  evolution.  In  future  work,  we  plan  to  in- 
vestigate the  evolution  of  the  packet  analytically  based 
on  Shrira’s  equations,  and  make  a comparison  with  our 
numerical  results. 
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Figure  13.  Contours  of  p/ po  from  a simulation  of  the  nonlinear  propagation  of  a wavepacket  with  initially  max 
dw/dz  « 0.761V.  Only  a portion  of  the  computational  frame  is  shown,  and  this  corresponds  to  a square  200/3  m 


on  each  side.  The  contour  increment  is  such  that  the  vertical  separation  between  unperturbed  isopycnals  is  2 m. 


The  times  corresponding  to  the  panels  are  (a)  2.5,  (b)  2.9  (c)  3.5,  (d)  3.9,  (e)  4.4  and  (f)  4.8  hr. 
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Conclusions 

We  have  attempted  with  some  very  idealized  numer- 
ical modelling  to  shed  some  light  on  the  structures  that 
may  be  responsible  for  the  variability  observed  in  the 
oceanic  thermocline.  This  attempt  has  been  successful 
to  a certain  extent.  We  were  able  to  show  a transition 
from  buoyancy  to  inertial  range  that  follows  the  fea- 
tures of  the  observed  spectra,  but  only  during  a break- 
ing event.  To  what  extent  this  simulated  transition  is 
related  to  the  actual  transition  in  the  oceanic  obser- 
vations is  still  not  clear.  We  were  able  to  give  three- 
dimensional  structure  to  a classical  shear  type  instabil- 
ity much  as  found  in  the  observations  and  at  about  the 
correct  vertical  scale.  In  addition,  we  described  an  un- 
expected kind  of  overturning  and  mixing  that  occurs  in 
regions  of  high  strain  rate.  The  correlation  for  overturns 
and  high  strain  rate  was  suggested  by  Alford  and  Pinkel 
(2000),  and  yet  it  is  not  clear  at  this  point  whether  the 
structures  we  were  able  to  simulate  are  directly  related 
to  any  of  the  overturns  reported  in  that  paper.  Fi- 
nally, two-dimensional  simulations  demonstrated  that 
a packet  with  approximately  the  correct  vertical  struc- 
ture can  propagate  through  a substantial  portion  of  a 
thermocline  as  observed  by  Alford  and  Pinkel  (2000) 
without  dispersing  radically  and  with  the  production  of 
overturns  on  the  scale  of  about  2 m which  is  entirely 
consistent  with  Alford  and  Pinkel’s  observations.  Un- 
fortunately, these  simulations  could  not  also  capture  the 
large  horizontal  scales  nor  the  rapid  advection  of  pack- 
ets that  is  suggested  by  Alford  and  Pinkel  (2000).  Thus, 
we  feel  that  there  has  been  some  progress  in  demon- 
strating the  possibility  of  analyzing  the  small-scale  fluc- 
tuations in  the  thermocline  with  numerical  simulations 
using  subgrid  scale  parameterizations,  and  look  forward 
to  advancing  this  work  so  that  closer  comparisons  with 
observations  may  be  possible. 
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